Global Sea Ice EDA and Analysis

Load in data

#load in global data set
globe_ice.df <- read.csv(file = "./data/global_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("global_ice" = "Ice",
         "global_anomaly" = "Anomaly")

summary(globe_ice.df)

##       Date           Date2                Year          Month      
##  Min.   :197901   Length:518         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200015                      Mean   :2000   Mean   : 6.48  
##  3rd Qu.:201105                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202202                      Max.   :2022   Max.   :12.00  
##    global_ice    global_anomaly 
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:   20   1st Qu.:   -1  
##  Median :   23   Median :    0  
##  Mean   :  -16   Mean   :  -39  
##  3rd Qu.:   25   3rd Qu.:    0  
##  Max.   :   28   Max.   :    2

#load in northern hemisphere data set
nohem_ice.df <- read.csv(file = "./data/nohem_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("nohem_ice" = "Value",
         "nohem_anomaly" = "Anomaly")

summary(nohem_ice.df)

##       Date           Date2                Year          Month      
##  Min.   :197901   Length:519         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200019                      Mean   :2000   Mean   : 6.47  
##  3rd Qu.:201106                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202203                      Max.   :2022   Max.   :12.00  
##    nohem_ice     nohem_anomaly  
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:    9   1st Qu.:   -1  
##  Median :   12   Median :    0  
##  Mean   :  -27   Mean   :  -39  
##  3rd Qu.:   14   3rd Qu.:    0  
##  Max.   :   16   Max.   :    1

#load in southern hemisphere data set
sohem_ice.df <- read.csv(file = "./data/sohem_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("sohem_ice" = "Value",
         "sohem_anomaly" = "Anomaly")

summary(sohem_ice.df)

##       Date           Date2                Year          Month      
##  Min.   :197901   Length:519         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200019                      Mean   :2000   Mean   : 6.47  
##  3rd Qu.:201106                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202203                      Max.   :2022   Max.   :12.00  
##    sohem_ice     sohem_anomaly  
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:    6   1st Qu.:    0  
##  Median :   12   Median :    0  
##  Mean   :  -27   Mean   :  -39  
##  3rd Qu.:   17   3rd Qu.:    0  
##  Max.   :   20   Max.   :    2

#load in hemispheric daily sea ice extent for duration calculations/plots
nohem_ice_daily.df <- read.csv(file = "./data/N_seaice_extent_daily_v3.0.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("nohem_extent" = "Extent") %>% 
  mutate(Date = as.Date(Date)) %>% 
  filter(!Date < "1988-12-31" & Year != 2022)  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data
  

summary(nohem_ice_daily.df)

##       Date                 Year          Month            Day      
##  Min.   :1988-12-31   Min.   :1988   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.:1997-04-01   1st Qu.:1997   1st Qu.: 4.00   1st Qu.: 8.0  
##  Median :2005-07-01   Median :2005   Median : 7.00   Median :16.0  
##  Mean   :2005-07-01   Mean   :2005   Mean   : 6.52   Mean   :15.7  
##  3rd Qu.:2013-09-30   3rd Qu.:2013   3rd Qu.:10.00   3rd Qu.:23.0  
##  Max.   :2021-12-31   Max.   :2021   Max.   :12.00   Max.   :31.0  
##   nohem_extent  
##  Min.   : 3.34  
##  1st Qu.: 8.30  
##  Median :11.86  
##  Mean   :11.11  
##  3rd Qu.:14.09  
##  Max.   :16.25

sohem_ice_daily.df <- read.csv(file = "./data/S_seaice_extent_daily_v3.0.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("sohem_extent" = "Extent") %>% 
  mutate(Date = as.Date(Date)) %>% 
  filter(!Date < "1988-12-31") %>%  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation
  filter(!Date < "1988-12-31" & Year != 2022)  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data

summary(sohem_ice_daily.df)

##       Date                 Year          Month            Day      
##  Min.   :1988-12-31   Min.   :1988   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.:1997-04-01   1st Qu.:1997   1st Qu.: 4.00   1st Qu.: 8.0  
##  Median :2005-07-01   Median :2005   Median : 7.00   Median :16.0  
##  Mean   :2005-07-01   Mean   :2005   Mean   : 6.52   Mean   :15.7  
##  3rd Qu.:2013-09-30   3rd Qu.:2013   3rd Qu.:10.00   3rd Qu.:23.0  
##  Max.   :2021-12-31   Max.   :2021   Max.   :12.00   Max.   :31.0  
##   sohem_extent  
##  Min.   : 2.08  
##  1st Qu.: 6.09  
##  Median :12.59  
##  Mean   :11.68  
##  3rd Qu.:17.24  
##  Max.   :20.20

#load in sea level data
sea_level.df <- read.csv(file = "./data/epa-sea-level.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("mean_sea_level" = "CSIRO.Adjusted.Sea.Level",
         "lowerCI_sea_level" = "Lower.Error.Bound",
         "upperCI_sea_level" = "Upper.Error.Bound")

summary(sea_level.df)

##      Date                Year      mean_sea_level lowerCI_sea_level
##  Length:35          Min.   :1979   Min.   :5.36   Min.   :5.10     
##  Class :character   1st Qu.:1988   1st Qu.:6.16   1st Qu.:5.91     
##  Mode  :character   Median :1996   Median :6.67   Median :6.39     
##                     Mean   :1996   Mean   :6.99   Mean   :6.72     
##                     3rd Qu.:2004   3rd Qu.:7.75   3rd Qu.:7.47     
##                     Max.   :2013   Max.   :9.33   Max.   :8.99     
##  upperCI_sea_level
##  Min.   :5.63     
##  1st Qu.:6.41     
##  Median :6.94     
##  Mean   :7.26     
##  3rd Qu.:8.03     
##  Max.   :9.66

#load in temp data
global_temp.df <- read.csv(file = "./data/monthly_globaltemp.csv", header = TRUE) %>% 
  as_tibble() %>% 
  arrange(-row_number()) %>%  #reverse row order so earliest years come first, matches sea ice df
  mutate(Date = as.Date(Date) %>% 
           as.yearmon()) #change to month year format

#join the three dataframes
total_ice.df <- inner_join(globe_ice.df, nohem_ice.df,
                           by = c("Date" = "Date",
                                  "Date2" = "Date2",
                                  "Year" = "Year",
                                  "Month" = "Month")) 

total_ice.df <- inner_join(total_ice.df, sohem_ice.df,
                           by = c("Date" = "Date",
                                  "Date2" = "Date2",
                                  "Year" = "Year",
                                  "Month" = "Month"))

#tidy up data
total_ice.df <- total_ice.df %>% 
  filter(global_ice != -9999) %>% #remove rows with -9999 (no value)
  dplyr::select(!Date) %>% 
  rename(Date = Date2) %>% 
  mutate(Date = str_replace_all(string = Date, 
                            pattern = "/", 
                            replacement = "-"))

total_ice.df <- total_ice.df %>% 
  mutate(Date = mdy(Date))
  

#check how summary looks now
summary(total_ice.df)

##       Date                 Year          Month         global_ice  
##  Min.   :1979-01-01   Min.   :1979   Min.   : 1.00   Min.   :16.3  
##  1st Qu.:1989-11-23   1st Qu.:1989   1st Qu.: 3.00   1st Qu.:20.3  
##  Median :2000-08-16   Median :2000   Median : 6.00   Median :23.5  
##  Mean   :2000-08-03   Mean   :2000   Mean   : 6.48   Mean   :23.0  
##  3rd Qu.:2011-05-08   3rd Qu.:2011   3rd Qu.: 9.00   3rd Qu.:25.3  
##  Max.   :2022-02-01   Max.   :2022   Max.   :12.00   Max.   :27.8  
##  global_anomaly     nohem_ice     nohem_anomaly      sohem_ice    
##  Min.   :-3.720   Min.   : 3.57   Min.   :-3.020   Min.   : 2.16  
##  1st Qu.:-0.725   1st Qu.: 8.56   1st Qu.:-0.843   1st Qu.: 5.82  
##  Median :-0.100   Median :12.10   Median :-0.185   Median :12.12  
##  Mean   :-0.279   Mean   :11.39   Mean   :-0.274   Mean   :11.57  
##  3rd Qu.: 0.390   3rd Qu.:14.30   3rd Qu.: 0.340   3rd Qu.:17.02  
##  Max.   : 1.650   Max.   :16.34   Max.   : 1.260   Max.   :19.76  
##  sohem_anomaly    
##  Min.   :-2.1300  
##  1st Qu.:-0.3200  
##  Median :-0.0350  
##  Mean   :-0.0059  
##  3rd Qu.: 0.3300  
##  Max.   : 1.8500

str(total_ice.df)

## tibble [516 × 9] (S3: tbl_df/tbl/data.frame)
##  $ Date          : Date[1:516], format: "1979-01-01" "1979-02-01" ...
##  $ Year          : int [1:516] 1979 1979 1979 1979 1979 1979 1979 1979 1979 1979 ...
##  $ Month         : int [1:516] 1 2 3 4 5 6 7 8 9 10 ...
##  $ global_ice    : num [1:516] 20.8 19.3 20.3 22.9 24.7 ...
##  $ global_anomaly: num [1:516] 1.39 0.95 0.88 1.4 1.24 1.6 1.4 0.82 0.35 0.13 ...
##  $ nohem_ice     : num [1:516] 15.4 16.2 16.3 15.4 13.9 ...
##  $ nohem_anomaly : num [1:516] 0.99 0.88 0.91 0.76 0.57 0.76 0.84 0.84 0.64 0.4 ...
##  $ sohem_ice     : num [1:516] 5.4 3.14 4 7.49 10.83 ...
##  $ sohem_anomaly : num [1:516] 0.4 0.07 -0.03 0.64 0.67 0.84 0.56 -0.02 -0.29 -0.28 ...

Plotting temporal trends in global and hemispheric sea ice change

total_ice_bar.df <- total_ice.df %>% 
  pivot_longer(cols = ends_with("ice"), #taking from wide to long
               names_to = "ice_type", #category = type of event 
               values_to = "area") %>%  #measurement = time till event
  filter(ice_type != "global_hem") %>% 
  mutate(ice_type = as.factor(ice_type) %>% 
           fct_recode("Northern" = "nohem_ice",
                      "Southern" = "sohem_ice"))
  
#make ggplotly  
sea.ice.gg <- total_ice_bar.df %>% 
  ggplot(aes(x = Date, y = area, fill = ice_type)) +
  geom_bar(stat = "identity", position = "stack", size = 1.5) +
  scale_y_continuous(limits = c(0,25)) +
  scale_x_date( #The date equiv of scale_x_continuous
    date_breaks = "5 year", #Breaks
    date_minor_breaks = "1 year", #Little lines in between
    date_labels = "%b '%y" #Label for the date,
  ) +
  xlab("") + ylab("Ice Area (millions km^2)") +
  scale_fill_manual(values = met.brewer("Hokusai3", direction = 1))
  
  

#ggplotly() + tooltip to control label
ggplotly(sea.ice.gg, tooltip = "text") #tooltip controls text hover/label

total_ice_line.df <- total_ice.df %>% 
  pivot_longer(cols = ends_with("ice"), #taking from wide to long
               names_to = "ice_type", #category = type of event 
               values_to = "extent") %>%  #measurement = time till event
  mutate(ice_type = as.factor(ice_type) %>% 
           fct_recode("Global" = "global_ice",
                      "Northern" = "nohem_ice",
                      "Southern" = "sohem_ice"))

Now make in plotly

#make plotly of sea ice
seaice.plotly <- plot_ly(total_ice_line.df, x = ~Date, y = ~extent, color = ~ice_type)
seaice.plotly <- seaice.plotly %>% add_lines()
seaice.plotly <- seaice.plotly %>% layout(
    xaxis = list(
      rangeselector = list(
        buttons = list(
          list(
            count = 6,
            label = "6 mo",
            step = "month",
            stepmode = "backward"),
          list(
            count = 1,
            label = "1 yr",
            step = "year",
            stepmode = "backward"),
          list(
          count = 5,
            label = "5 yr",
            step = "year",
            stepmode = "backward"),
          list(
          count = 10,
            label = "10 yr",
            step = "year",
            stepmode = "backward"),
          list(
            count = 1,
            label = "YTD",
            step = "year",
            stepmode = "todate"),
          list(step = "all"))),

      rangeslider = list(type = "date")),
    
    yaxis = list(title = "Sea Ice Extent (millions km^2)"))

seaice.plotly

Let’s make some plots looking at annual, as opposed to monthly, shifts in sea ice

total_ice_year.df <- total_ice_line.df %>% 
  group_by(Year, ice_type) %>% 
  summarize(
    mean_ice_extent = mean(extent)
  ) %>% 
  ungroup() %>% 
  filter(Year != 2022) #remove last year since incomplete, throws off annual mean

#now plot
seaice.year.plotly <- plot_ly(data = total_ice_year.df, x = ~Year, y = ~mean_ice_extent, color = ~ice_type)
seaice.year.plotly <- seaice.year.plotly %>% add_lines()
seaice.year.plotly <- seaice.year.plotly %>% layout(
  xaxis = list("Year"),
  yaxis = list(title = "Sea Ice Extent (millions km^2)"))

seaice.year.plotly

Not surprisingly, it looks like we have a declining trend in sea ice extent. Let’s create some simple linear models to compare trends.

#make global dataset
global_ice_year.df <- total_ice_year.df %>% 
  filter(ice_type %in% "Global")

#make northern dataset
nohem_ice_year.df <- total_ice_year.df %>% 
  filter(ice_type %in% "Northern")

#make southern dataset
sohem_ice_year.df <- total_ice_year.df %>% 
  filter(ice_type %in% "Southern")

#construct three simple linear models comparing year with ice extent to get annual mean trend
mod1 <- lm(mean_ice_extent ~ Year, data = global_ice_year.df)
mod2 <- lm(mean_ice_extent ~ Year, data = nohem_ice_year.df)
mod3 <- lm(mean_ice_extent ~ Year, data = sohem_ice_year.df)

#check out model summaries
summary(mod1)

## 
## Call:
## lm(formula = mean_ice_extent ~ Year, data = global_ice_year.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0262 -0.3830 -0.0058  0.3960  1.2481 
## 
## Coefficients:
##              Estimate Std. Error t value        Pr(>|t|)    
## (Intercept) 121.21108   12.45421    9.73 0.0000000000032 ***
## Year         -0.04912    0.00623   -7.89 0.0000000009479 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.507 on 41 degrees of freedom
## Multiple R-squared:  0.603,  Adjusted R-squared:  0.593 
## F-statistic: 62.2 on 1 and 41 DF,  p-value: 0.000000000948

summary(mod2)

## 
## Call:
## lm(formula = mean_ice_extent ~ Year, data = nohem_ice_year.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5010 -0.1687  0.0225  0.1707  0.3300 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 119.55461    5.33322    22.4 <0.0000000000000002 ***
## Year         -0.05409    0.00267   -20.3 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.217 on 41 degrees of freedom
## Multiple R-squared:  0.909,  Adjusted R-squared:  0.907 
## F-statistic:  411 on 1 and 41 DF,  p-value: <0.0000000000000002

summary(mod3)

## 
## Call:
## lm(formula = mean_ice_extent ~ Year, data = sohem_ice_year.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9825 -0.2006  0.0013  0.1675  1.0565 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.67151   10.20383    0.16     0.87
## Year         0.00496    0.00510    0.97     0.34
## 
## Residual standard error: 0.415 on 41 degrees of freedom
## Multiple R-squared:  0.0226, Adjusted R-squared:  -0.00127 
## F-statistic: 0.947 on 1 and 41 DF,  p-value: 0.336

When will sea ice disappear from the earth globally and in each hemisphere?

Based on these simple linear models, if mean annual sea ice trends continued at the same rate as they have since 1979, we would predict that average sea ice extent will reach 0 in the northern hemisphere by the year 2210, and it will reach 0 globally and, thus, in the southern hemisphere, by the year 2468.

Let’s calculate the duration of sea ice expansion and and recession for each year

#join both daily data frames
bothhem_ice_daily.df <- inner_join(nohem_ice_daily.df, sohem_ice_daily.df,
                           by = c("Date" = "Date",
                                  "Year" = "Year",
                                  "Month" = "Month",
                                  "Day" = "Day"))

#for northern hemisphere
bothhem_ice_daily.df <- bothhem_ice_daily.df %>% 
  mutate(nohem_ice_diff = nohem_extent - lag(nohem_extent),
         sohem_ice_diff = sohem_extent - lag(sohem_extent),
         nohem_ice_trend = case_when(
           nohem_ice_diff > 0 ~ "increasing",
           nohem_ice_diff %in% 0 ~ "no change",
           nohem_ice_diff < 0 ~ "decreasing"),
         sohem_ice_trend = case_when(
           sohem_ice_diff > 0 ~ "increasing",
           sohem_ice_diff %in% 0 ~ "no change",
           sohem_ice_diff < 0 ~ "decreasing")
  ) %>% 
  na.omit()
  
  
summary(bothhem_ice_daily.df)

##       Date                 Year          Month            Day      
##  Min.   :1989-01-01   Min.   :1989   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.:1997-04-02   1st Qu.:1997   1st Qu.: 4.00   1st Qu.: 8.0  
##  Median :2005-07-02   Median :2005   Median : 7.00   Median :16.0  
##  Mean   :2005-07-02   Mean   :2005   Mean   : 6.52   Mean   :15.7  
##  3rd Qu.:2013-10-01   3rd Qu.:2013   3rd Qu.:10.00   3rd Qu.:23.0  
##  Max.   :2021-12-31   Max.   :2021   Max.   :12.00   Max.   :31.0  
##   nohem_extent    sohem_extent   nohem_ice_diff    sohem_ice_diff   
##  Min.   : 3.34   Min.   : 2.08   Min.   :-0.3970   Min.   :-0.5180  
##  1st Qu.: 8.30   1st Qu.: 6.09   1st Qu.:-0.0570   1st Qu.:-0.0730  
##  Median :11.86   Median :12.59   Median :-0.0050   Median : 0.0200  
##  Mean   :11.11   Mean   :11.68   Mean   :-0.0001   Mean   :-0.0001  
##  3rd Qu.:14.09   3rd Qu.:17.24   3rd Qu.: 0.0540   3rd Qu.: 0.0880  
##  Max.   :16.25   Max.   :20.20   Max.   : 0.3630   Max.   : 0.3310  
##  nohem_ice_trend    sohem_ice_trend   
##  Length:12053       Length:12053      
##  Class :character   Class :character  
##  Mode  :character   Mode  :character  
##                                       
##                                       
## 

#now summarize number of increasing and decreasing days for each year
#northern hemisphere
nohem_ice_duration.df <- bothhem_ice_daily.df %>% 
  group_by(Year, nohem_ice_trend) %>% 
  summarise(
    no_days = length(nohem_ice_trend)
  ) %>% 
  mutate(nohem_ice_trend = as.factor(nohem_ice_trend)) %>% 
  ungroup()


#southern hemisphere
sohem_ice_duration.df <- bothhem_ice_daily.df %>% 
  group_by(Year, sohem_ice_trend) %>% 
  summarise(
    no_days = length(sohem_ice_trend) 
  ) %>% 
  mutate(sohem_ice_trend = as.factor(sohem_ice_trend)) %>% 
  ungroup()

Now make plots showing trends in days if increasing and decreasing ice

#northern hemisphere plot
nohem.ice.dur.plotly <- plot_ly(data = nohem_ice_duration.df, x = ~Year, y = ~no_days, color = ~nohem_ice_trend)
nohem.ice.dur.plotly <- nohem.ice.dur.plotly %>% add_lines()
nohem.ice.dur.plotly <- nohem.ice.dur.plotly %>% layout(
  xaxis = list("Year"),
  yaxis = list(title = "No. of Days"))

nohem.ice.dur.plotly

#southern hemisphere plot
sohem.ice.dur.plotly <- plot_ly(data = sohem_ice_duration.df, x = ~Year, y = ~no_days, color = ~sohem_ice_trend)
sohem.ice.dur.plotly <- sohem.ice.dur.plotly %>% add_lines()
sohem.ice.dur.plotly <- sohem.ice.dur.plotly %>% layout(
  xaxis = list("Year"),
  yaxis = list(title = "No. of Days"))

sohem.ice.dur.plotly

Now do the same but for anomaly data

#alter df for anomaly plot
total_ice_anomaly.df <- total_ice.df %>% 
  pivot_longer(cols = ends_with("anomaly"), #taking from wide to long
               names_to = "ice_type", #category = type of event 
               values_to = "anomaly") %>%  #measurement = time till event
  mutate(ice_type = as.factor(ice_type) %>% 
           fct_recode("Global" = "global_anomaly",
                      "Northern" = "nohem_anomaly",
                      "Southern" = "sohem_anomaly"))

#make plotly of sea ice anomaly values
iceanomaly.plotly <- plot_ly(total_ice_anomaly.df, x = ~Date, y = ~anomaly, color = ~ice_type)
iceanomaly.plotly <- iceanomaly.plotly %>% add_lines()
iceanomaly.plotly <- iceanomaly.plotly %>% layout(
    xaxis = list(
      rangeselector = list(
        buttons = list(
          list(
            count = 6,
            label = "6 mo",
            step = "month",
            stepmode = "backward"),
          list(
            count = 1,
            label = "1 yr",
            step = "year",
            stepmode = "backward"),
          list(
          count = 5,
            label = "5 yr",
            step = "year",
            stepmode = "backward"),
          list(
          count = 10,
            label = "10 yr",
            step = "year",
            stepmode = "backward"),
          list(
            count = 1,
            label = "YTD",
            step = "year",
            stepmode = "todate"),
          list(step = "all"))),

      rangeslider = list(type = "date")),
    
    yaxis = list(title = "Sea Ice Anomalies"))

iceanomaly.plotly

It looks like global anomalies in sea ice are mostly driven by the northern hemisphere sea ice as opposed to the southern hemisphere. Let’s rtun some simple linear models to check it out.

#global anomaly vs southern hemisphere
mod4 <- lm(global_anomaly ~ sohem_anomaly, data = total_ice.df)
summary(mod4)

## 
## Call:
## lm(formula = global_anomaly ~ sohem_anomaly, data = total_ice.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.692 -0.546  0.113  0.594  1.571 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)    -0.2739     0.0348   -7.87    0.000000000000021 ***
## sohem_anomaly   0.8863     0.0625   14.19 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.791 on 514 degrees of freedom
## Multiple R-squared:  0.281,  Adjusted R-squared:  0.28 
## F-statistic:  201 on 1 and 514 DF,  p-value: <0.0000000000000002

#global anomaly vs northern hemisphere
mod5 <- lm(global_anomaly ~ nohem_anomaly, data = total_ice.df)
summary(mod5)

## 
## Call:
## lm(formula = global_anomaly ~ nohem_anomaly, data = total_ice.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1877 -0.3161 -0.0229  0.3070  1.8240 
## 
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)    -0.0205     0.0259   -0.79                0.43    
## nohem_anomaly   0.9434     0.0309   30.49 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.557 on 514 degrees of freedom
## Multiple R-squared:  0.644,  Adjusted R-squared:  0.643 
## F-statistic:  930 on 1 and 514 DF,  p-value: <0.0000000000000002

#plots of both models above
#sohem
total_ice.df %>% 
  ggplot(aes(x = sohem_anomaly,  y = global_anomaly)) +
  stat_smooth(formula = "y ~ x",
              method = "lm")

#nohem
total_ice.df %>% 
  ggplot(aes(x = nohem_anomaly,  y = global_anomaly)) +
  stat_smooth(formula = "y ~ x",
              method = "lm")

As expected from an earlier visual assessment of anomaly data, the model with northern hemisphere anomaly explains 64.3% of variation in global anomaly data, relative to only 28% explained by the southern hemisphere anomalies.

Modeling relationship between sea ice change and sea level rise

#prep sea ice df for comparison with sea level data
total_ice_march.df <- total_ice.df %>% 
  filter(Month == 3) %>% #only use ice data from March, matches collectiion of annual sea level data
  dplyr::select(Year, global_ice, global_anomaly) %>% 
  mutate(global_ice_change = global_ice - lag(global_ice)) %>% #create annual change variable
  filter(Year != 1979) #remove first year with missing change value

#remove Date column
sea_level.df <- sea_level.df %>% 
  dplyr::select(!Date)

#join dataframes
sea_levelandice.df <- inner_join(total_ice_march.df, sea_level.df,
                                 by = ("Year" = "Year"))

#check out summary
summary(sea_levelandice.df)

##       Year        global_ice   global_anomaly    global_ice_change
##  Min.   :1980   Min.   :17.6   Min.   :-1.8300   Min.   :-1.1500  
##  1st Qu.:1988   1st Qu.:19.1   1st Qu.:-0.3375   1st Qu.:-0.5875  
##  Median :1996   Median :19.4   Median :-0.0100   Median :-0.1950  
##  Mean   :1996   Mean   :19.4   Mean   :-0.0253   Mean   :-0.0085  
##  3rd Qu.:2005   3rd Qu.:19.9   3rd Qu.: 0.4525   3rd Qu.: 0.4750  
##  Max.   :2013   Max.   :20.5   Max.   : 1.0000   Max.   : 2.0900  
##  mean_sea_level lowerCI_sea_level upperCI_sea_level
##  Min.   :5.60   Min.   :5.34      Min.   :5.85     
##  1st Qu.:6.17   1st Qu.:5.92      1st Qu.:6.42     
##  Median :6.73   Median :6.46      Median :7.00     
##  Mean   :7.04   Mean   :6.76      Mean   :7.30     
##  3rd Qu.:7.75   3rd Qu.:7.48      3rd Qu.:8.04     
##  Max.   :9.33   Max.   :8.99      Max.   :9.66

#manual colors for legend
colors_icelevel <- c("Global Sea Level" = "#3399FF", "Global Sea Ice" = "#993399")

#plot sea ice and level together
sea_levelandice.df %>% 
  ggplot(aes(x = Year)) +
  geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
  geom_line(aes(y = mean_sea_level*3, color = "Global Sea Level")) +
  geom_ribbon(aes(ymin = lowerCI_sea_level*3, ymax = upperCI_sea_level*3), fill = "#99CCFF", alpha = 0.25) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Extent (millions km^2)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./3, name = "Mean Global Sea Level Change (in)")) +
  scale_color_manual("", #manual legend
                     values = colors_icelevel) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))

#construct simple linear regression models
mod6 <- lm(mean_sea_level ~ global_ice, data = sea_levelandice.df)
mod7 <- lm(mean_sea_level ~ global_ice_change + poly(global_ice_change, 2), data = sea_levelandice.df)

#model summary
summary(mod6)

## 
## Call:
## lm(formula = mean_sea_level ~ global_ice, data = sea_levelandice.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.484 -0.823 -0.221  0.589  2.437 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   15.868      5.681    2.79   0.0087 **
## global_ice    -0.454      0.292   -1.56   0.1297   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.07 on 32 degrees of freedom
## Multiple R-squared:  0.0703, Adjusted R-squared:  0.0412 
## F-statistic: 2.42 on 1 and 32 DF,  p-value: 0.13

summary(mod7)

## 
## Call:
## lm(formula = mean_sea_level ~ global_ice_change + poly(global_ice_change, 
##     2), data = sea_levelandice.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.859 -0.656 -0.208  0.559  2.352 
## 
## Coefficients: (1 not defined because of singularities)
##                             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)                    7.037      0.170   41.31 <0.0000000000000002 ***
## global_ice_change              0.116      0.210    0.55              0.5839    
## poly(global_ice_change, 2)1       NA         NA      NA                  NA    
## poly(global_ice_change, 2)2    2.866      0.993    2.89              0.0071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.993 on 31 degrees of freedom
## Multiple R-squared:  0.218,  Adjusted R-squared:  0.167 
## F-statistic: 4.32 on 2 and 31 DF,  p-value: 0.0222

#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod6)

par(mfrow = c(2,2))
plot(mod7)

#plot relationships
sea_levelandice.df %>% 
  ggplot(aes(x = global_ice,  y = mean_sea_level)) +
  stat_smooth(formula = "y ~ poly(x,2)",
              method = "lm")

sea_levelandice.df %>% 
  ggplot(aes(x = global_ice_change,  y = mean_sea_level)) +
  stat_smooth(formula = "y ~ poly(x,2)",
              method = "lm")

Monthly sea ice and temperature

#convert date to month year to match temperature data
total_ice_yearmon.df <- total_ice.df %>% 
  mutate(Date = as.yearmon(Date))

#split  temp df into two dfsbased on source
global_temp_gis.df <- global_temp.df %>% 
  filter(Source %in% "GISTEMP")

global_temp_gcag.df <- global_temp.df %>% 
  filter(Source %in% "GCAG")

#join dataframes
ice_temp_gis.df <- inner_join(total_ice_yearmon.df, global_temp_gis.df,
                          by = ("Date" = "Date")) 

ice_temp_gcag.df <- inner_join(total_ice_yearmon.df, global_temp_gcag.df,
                          by = ("Date" = "Date"))

#manual colors for legend
colors_icetemp <- c("Global Temperature" = "#3399FF", "Global Sea Ice" = "#993399")

#plot sea ice and level together
ice_temp_gis.df %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
  geom_line(aes(y = Mean*25, color = "Global Temperature")) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Extent (millions km^2)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./25, name = "Mean Global Temp Change (degrees Celsius)")) +
  scale_color_manual("", #manual legend
                     values = colors_icetemp) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))

 #plot sea ice and level together
ice_temp_gis.df %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = global_anomaly, color = "Global Sea Ice")) +
  geom_line(aes(y = Mean*3, color = "Global Temperature")) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Anomalies",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./3, name = "Mean Global Temp Change (degrees Celsius)")) +
  scale_color_manual("", #manual legend
                     values = colors_icetemp) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))                   

#construct simple linear regression models
mod8 <- lm(global_ice ~ Mean, data = ice_temp_gis.df)
mod9 <- lm(global_anomaly ~ Mean, data = ice_temp_gis.df)

#model summary
summary(mod8)

## 
## Call:
## lm(formula = global_ice ~ Mean, data = ice_temp_gis.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.227 -2.568  0.848  2.146  4.510 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   24.383      0.287   85.01 < 0.0000000000000002 ***
## Mean          -2.571      0.546   -4.71            0.0000033 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.75 on 452 degrees of freedom
## Multiple R-squared:  0.0468, Adjusted R-squared:  0.0447 
## F-statistic: 22.2 on 1 and 452 DF,  p-value: 0.00000327

summary(mod9)

## 
## Call:
## lm(formula = global_anomaly ~ Mean, data = ice_temp_gis.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8921 -0.3729  0.0605  0.4468  1.6615 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   0.6942     0.0694    10.0 <0.0000000000000002 ***
## Mean         -1.6366     0.1319   -12.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.666 on 452 degrees of freedom
## Multiple R-squared:  0.254,  Adjusted R-squared:  0.252 
## F-statistic:  154 on 1 and 452 DF,  p-value: <0.0000000000000002

#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod8)

par(mfrow = c(2,2))
plot(mod9)

#plot relationships
ice_temp_gis.df %>% 
  ggplot(aes(x = Mean,  y = global_ice)) +
  stat_smooth(formula = "y ~ x",
              method = "lm")

ice_temp_gis.df %>% 
  ggplot(aes(x = Mean,  y = global_anomaly)) +
  stat_smooth(formula = "y ~ x", 
              method = "lm")

Spatial plots of sea ice

Read in a few csvs to get started

#load shapefiles
seaice.jan1980 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp")

## Reading layer `extent_N_198001_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 143 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3400000 ymin: -4875000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North

seaice.jan2020 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202001_polygon_v3.0.shp")

## Reading layer `extent_N_202001_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202001_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3375000 ymin: -4925000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North

seaice.sep1980 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp")

## Reading layer `extent_N_198009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 118 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2550000 ymin: -3250000 xmax: 2350000 ymax: 2100000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North

seaice.sep2020 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp")

## Reading layer `extent_N_202009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2525000 ymin: -2800000 xmax: 2350000 ymax: 4950000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North

# seaice.jan1980.transform <- stTransform(seaice.jan1980, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

#make basic plot of layers
ggplot() +
  geom_sf(data = seaice.jan1980, fill = "blue") +
  geom_sf(data = seaice.jan2020, fill = "purple")

ggplot() +
  geom_sf(data = seaice.sep1980, fill = "red") +
  geom_sf(data = seaice.sep2020, fill = "green")

ggplot() +
  geom_sf(data = seaice.jan2020, fill = "brown") +
  geom_sf(data = seaice.sep2020, fill = "yellow")

#load shapefile
seaice.jan1980 <- readOGR("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp", layer: "extent_N_198001_polygon_v3.0"
## with 143 features
## It has 1 fields
## Integer64 fields read as strings:  FID

#check class and crs
class(seaice.jan1980)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

crs(seaice.jan1980)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +a=6378273
## +rf=298.279411123064 +units=m +no_defs +type=crs 
## WKT2 2019 representation:
## PROJCRS["NSIDC Sea Ice Polar Stereographic North",
##     BASEGEOGCRS["Unspecified datum based upon the Hughes 1980 ellipsoid",
##         DATUM["Not specified (based on Hughes 1980 ellipsoid)",
##             ELLIPSOID["Hughes 1980",6378273,298.279411123064,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4054]],
##     CONVERSION["US NSIDC Sea Ice polar stereographic north",
##         METHOD["Polar Stereographic (variant B)",
##             ID["EPSG",9829]],
##         PARAMETER["Latitude of standard parallel",70,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8832]],
##         PARAMETER["Longitude of origin",-45,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8833]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",south,
##             MERIDIAN[45,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",south,
##             MERIDIAN[135,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Polar research."],
##         AREA["Northern hemisphere - north of 60°N onshore and offshore, including Arctic."],
##         BBOX[60,-180,90,180]],
##     ID["EPSG",3411]]

extent(seaice.jan1980)

## class      : Extent 
## xmin       : -3400000 
## xmax       : 3600000 
## ymin       : -4875000 
## ymax       : 5850000

seaice.jan1980.WGS84 <- spTransform(seaice.jan1980, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

crs(seaice.jan1980.WGS84)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

extent(seaice.jan1980.WGS84)

## class      : Extent 
## xmin       : -180 
## xmax       : 179.7 
## ymin       : 38.01 
## ymax       : 83.67

broom::tidy(seaice.jan1980.WGS84)

## Regions defined for each Polygons

## # A tibble: 4,470 × 7
##     long   lat order hole  piece group id   
##    <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
##  1  125.  38.8     1 FALSE 1     0.1   0    
##  2  125.  38.7     2 FALSE 1     0.1   0    
##  3  124.  39.4     3 FALSE 1     0.1   0    
##  4  125.  39.5     4 FALSE 1     0.1   0    
##  5  125.  39.5     5 FALSE 1     0.1   0    
##  6  125.  39.7     6 FALSE 1     0.1   0    
##  7  125.  39.7     7 FALSE 1     0.1   0    
##  8  125.  39.5     8 FALSE 1     0.1   0    
##  9  125.  39.6     9 FALSE 1     0.1   0    
## 10  125.  39.2    10 FALSE 1     0.1   0    
## # … with 4,460 more rows

Make basemap

#manual bounding box
seaice.box <- matrix(c(-180, 180, 38.01, 83.67),
                     ncol  = 2,
                     byrow = TRUE)
colnames(seaice.box) <- c("min", "max")
rownames(seaice.box) <- c("x", "y")
seaice.box

##       min    max
## x -180.00 180.00
## y   38.01  83.67

#Get the base map (foundational layer)
seaice_base.map <- get_map(
                       location = seaice.box,
                       source   = "stamen",
                       maptype  = "watercolor",
                       crop = TRUE
                      )

## Source : http://tile.stamen.com/terrain/3/0/0.png

## Source : http://tile.stamen.com/terrain/3/1/0.png

## Source : http://tile.stamen.com/terrain/3/2/0.png

## Source : http://tile.stamen.com/terrain/3/3/0.png

## Source : http://tile.stamen.com/terrain/3/4/0.png

## Source : http://tile.stamen.com/terrain/3/5/0.png

## Source : http://tile.stamen.com/terrain/3/6/0.png

## Source : http://tile.stamen.com/terrain/3/7/0.png

## Source : http://tile.stamen.com/terrain/3/8/0.png

## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/0.png.
## Source : http://tile.stamen.com/terrain/3/0/1.png
## Source : http://tile.stamen.com/terrain/3/1/1.png
## Source : http://tile.stamen.com/terrain/3/2/1.png
## Source : http://tile.stamen.com/terrain/3/3/1.png
## Source : http://tile.stamen.com/terrain/3/4/1.png
## Source : http://tile.stamen.com/terrain/3/5/1.png
## Source : http://tile.stamen.com/terrain/3/6/1.png
## Source : http://tile.stamen.com/terrain/3/7/1.png
## Source : http://tile.stamen.com/terrain/3/8/1.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/1.png.
## Source : http://tile.stamen.com/terrain/3/0/2.png
## Source : http://tile.stamen.com/terrain/3/1/2.png
## Source : http://tile.stamen.com/terrain/3/2/2.png
## Source : http://tile.stamen.com/terrain/3/3/2.png
## Source : http://tile.stamen.com/terrain/3/4/2.png
## Source : http://tile.stamen.com/terrain/3/5/2.png
## Source : http://tile.stamen.com/terrain/3/6/2.png
## Source : http://tile.stamen.com/terrain/3/7/2.png
## Source : http://tile.stamen.com/terrain/3/8/2.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/2.png.
## Source : http://tile.stamen.com/terrain/3/0/3.png
## Source : http://tile.stamen.com/terrain/3/1/3.png
## Source : http://tile.stamen.com/terrain/3/2/3.png
## Source : http://tile.stamen.com/terrain/3/3/3.png
## Source : http://tile.stamen.com/terrain/3/4/3.png
## Source : http://tile.stamen.com/terrain/3/5/3.png
## Source : http://tile.stamen.com/terrain/3/6/3.png
## Source : http://tile.stamen.com/terrain/3/7/3.png
## Source : http://tile.stamen.com/terrain/3/8/3.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/3.png.

ggmap(seaice_base.map) +
  geom_polygon(data = seaice.jan1980.WGS84, aes(x = long, y = lat))

## Regions defined for each Polygons

leaflet(data = seaice.jan1980.WGS84) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons()